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In a recent series of scanning probe experiments, it became possible to visualize local electron 
flow in a two-dimensional electron gas. In this paper, a Green's function technique is presented 
that enables efficient calculation of the quantity measured in such experiments. Efficient means 
that the computational effort scales like M 3 N (M is the width of the tight-binding lattice used, 
and N is its length), which is a factor MN better than the standard recursive technique for the 
same problem. Moreover, within our numerical framework it is also possible to calculate (with 
the same computational effort M 3 N) the local density of states, the electron density, and the 
current distribution in the sample, which are not accessible with the standard recursive method. 
Furthermore, an imaging method is discussed where the scanning tip can be used to measure the 
local chemical potential. The numerical technique is used to study electron flow through a quantum 
point contact. All features seen in experiments on this system are reproduced and a new interference 
effect is observed resulting from the crossing of coherent beams of electron flow. 



I. INTRODUCTION 

Electronic transport properties in mesoscopic systems 
have gained a lot of attention in the last two decades 
(for a review see e.g. Refs. 00). A lot of effort al- 
ready went into the study of global transport quantities 
(mainly conductance), and a range of interesting phe- 
nomena emerged: e.g. quantized conductance in a point 
contact 3 , quantum Hall effect 4 , and universal conduc- 
tance fluctuations^ to name a few. 

Recently, scanning probe methods have offered the in- 
teresting possibility of obtaining local information on 
electron transport, not accessible with ordinary global 
measurements. For example, by using the tip of an STM 
as a localized scatterer of electrons and measuring the 
conductance change of the sample as a function of tip 
position, it has been possible to obtain a spatial map of 
electron flow in a two-dimensional electron gas§i&. Such 
visualizations of current flow are interesting both from 
an experimental and a theoretical point of view. 

The results obtained in these experiments were mainly 
interpreted with the help of electron or flux density cal- 
culations. Here, a numerical technique is developed that 
allows to calculate in an efficient way the quantity that 
is actually measured, which is a conductance difference 
as a function of the tip position. The standard recur- 
sive Green's function method used in Ref. to obtain 
the same result is not very efficient because one has to 
start over the complete conductance calculation for ev- 
ery position of the tip. The numerical effort then scales 
like M 4 N 2 in the number of operations, where M is the 
width and N is the length of the system (N ^ 1), and 
one therefore is limited to small systems. On the other 
hand, our method scales like M 3 N for the same prob- 
lem. Moreover, with our method it is straightforward to 
obtain other relevant information like the local density 
of states and the electron and current density distribu- 
tions which can then be compared with the experimental 
observables. All these quantities can be calculated with 
the same computational effort scaling like M 3 N, whereas 



they are unavailable with the standard technique. 

The imaging technique with the tip as a scatterer will 
not always give the expected results in the regime of high 
magnetic fields. To obtain more information about elec- 
tron flow in this high field range, a measuring method is 
proposed where the STM tip is used as a local voltage 
probe. This method can also be treated within the pro- 
posed numerical framework, with a computational effort 
scaling again as M 3 N. 

The paper is subdivided as follows, fn the next section, 
the imaging method used in Refs. |fJl3 is briefly reviewed 
and the voltage probe method is described. In section lTTTl 
we discuss how the considered systems can be described 
within a tight-binding model. Subsequently, the numeri- 
cal method is presented in some detail in section Hvl and 
in the appendices. Finally, the method is applied to study 
electron flow through a quantum point contact. 



II. IMAGING METHODS 
A. Tip as a Local Scatterer 

The experiments in Refs. 00 take place as follows: cur- 
rent is passed via two Ohmic contacts through a two- 
dimensional electron gas, while simultaneously a scan- 
ning probe tip is moved across the sample. The local elec- 
trostatic potential that results from a negative voltage on 
the tip can function as a scattering center for electrons in 
the device. If the tip is placed over a region where a lot 
of electrons are flowing, the conductance of the sample 
will decrease considerably because of enhanced backscat- 
tering due to the tip, whereas in an area of low electron 
flow, the conductance decrease will be small. As such, by 
mapping the conductance decrease of the sample to the 
tip position, one gets a spatial profile of electron flow. 

The relevant quantity for this imaging method is thus 
the position dependent conductance difference: 



Ag(x,y) = g -g t (x,y), 



(1) 
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where gt(x, y) is the conductance with the tip positioned 
on a point (x, y) and go is the conductance without 
the tip. Conductances can be calculated within the 
Landauer-Biittiker formalism, which relates conductance 
to the electron transmission coefficient T between the 
contacts on opposite sides of the sample: 

2e 2 

Please note that since g is a two-terminal conductance, it 
is symmetric with respect to reversal of the direction of an 
applied magnetic field: </(+B) = g(— B) 10 , and therefore 
the map of electron flow obtained with this scattering 
method has the same symmetry. 

B. Tip as a Voltage Probe 

The imaging method in the previous section can give 
very nice visualizations of the electron flow, as proven by 
the experimental results of electron transport through a 
quantum point contact^. One of the interesting subjects 
to study with this imaging method would be the quantum 
Hall effect, because even after intensive theoretical and 
experimental effort since its discovery^, some of the de- 
tails of electron transport in this regime remain unclear. 
Unfortunately, the appearance of edge states in the high 
field regime leads to the suppression of backscattering 11 , 
so that the method described above will not yield the 
expected results in the quantum Hall limit. 

However, a picture of this edge-state transport can be 
obtained by investigating the local chemical potential in 
the sample. This can be done using the STM tip as a 
voltage probe; the STM tip voltage equilibrates itself to 
the local chemical potential when electrons are allowed 
to tunnel into and out of the tip. Experimentally, this 
technique has already been used to probe the potential 
distribution at metal-insulator-metal interfaces^ and at 
grain boundaries 13 . Gramespacher et al. have studied 
this measurement mode theoretically, and found Green's 
function expressions for the sample-tip transmission co- 
efficients which generalize the Bardeen formula^. For 
weak coupling between the probe and the sample, these 
transmission coefficients can be expressed in terms of so- 
called partial densities of states^. 

In the following, we will however start again from the 
standard Landauer-Biittiker formalism, but now one has 
three contacts: the left and right contact through which 



a fixed current is passed, and the STM tip that measures 
a voltage. For this 3-lead structure at temperature T = 
0, one obtains the following expression for the voltage 
measured on the tipi: 

Kip -V L = - Tt '^ L (V R V L ), (3) 

with Vp the voltage on lead p. Our numerical method 
thus should be able to calculate the transmission coef- 
ficients Itip^rm from the left and right contact to the 
STM tip, as a function of the position of the tip. 

It should be indicated that this imaging mode corre- 
sponds to making a three-terminal measurement. This 
means that the spatial map one obtains will not be in- 
variant under the reversal of the direction of an applied 
magnetic field. The only symmetry relations that hold 
are: T pq (+B) = T qp (-B) 10 . Therefore, it should be clear 
that the information obtained when the tip is used as a 
voltage probe is different from that when the tip is used 
as a local scatterer, and as such both imaging methods 
will contribute differently to our understanding of elec- 
tron flow in 2D systems. 



C. Charge and Current Density 

One can (at least in principal) experimentally observe 
the total current density in the sample, by measuring the 
magnetic field it induces. Another physically observable 
quantity is the total charge density, due to the electro- 
static field it generates. Both of these quantities can also 
be calculated within the numerical framework presented 
in the next sections. 

Since we have access to the electron density, it would 
also be possible to include self-consistently the effect of 
the applied voltage and the resultant current flow on the 
transport properties. However, this will not be pursued 
further in the current paper. 



III. SYSTEM MODELLING 

The system we will consider is a two-dimensional de- 
vice connected to two semi-infinite leads extending in the 
x-direction. By discretizing the Schrodinger equation, 
one obtains the standard nearest neighbor tight-binding 
model for the system Hamiltonian: 



oo M oo M—X 

H = ^2 ^2 £mn\m,n)(m,n\ + ^ ^ (t^ Jm,n + l)(m,n\ + t^ Jm + l,n) (m,n\ + h.c.) . (4) 
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where (m, n) label the sites on the lattice, and e mn are the resent atoms as in the usual tight-binding model; rather 
on-site energies. Note that the sites are not meant to rep- 
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they may represent a region containing many atoms, but 
this region should be small with respect to physically 
relevant quantities such as the Fermi wavelength. The 
quantities t^, and t v m n in Eq. Q give the hopping am- 
plitude in the horizontal, respectively vertical direction. 
A square lattice is assumed, so in the absence of a mag- 
netic field one has: 



t x =t v = -t = - 

b m,n m,n 2m * q2 : 



(5) 



with a the lattice spacing, and m* the effective mass of 
the electron. 

The leads extend from n — — 00..O for the left lead, 
and n = N + 1, .., 00 for the right one. The device itself 
is comprised by the other M x N sites. In the leads, 
only homogeneous fields perpendicular to the 2D device 
will be considered, while both homogeneous and inho- 
mogeneous fields can be treated in the central device. 
The fields are described by Peierls substitution 1 ?, which 
changes the hopping parameters as: 



x{y) 



_ t e -ie/hJ A.dl 



(6) 



where J A . dl is the integral of the vector potential along 
the hopping path. 



IV. NUMERICAL METHOD 
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FIG. 1: View on the Green's functions necessary to calculate 
the relevant physical quantities in our problem. The first and 
last column of the tight binding lattice are colored grey to 
indicate that they are influenced by the self-energies El and 
E_r of the leads. 



of operations that scales like M 3 N. As for the self- 
energies of the leads, they can be calculated in sev- 
eral ways: in the absence of a magnetic field they are 
known analytically^, while in the case of a homogeneous 
magnetic field one needs to resort to numerical meth- 
ods, for instance based on the solution of the so-called 
Harper's equation^. However, we have chosen for an- 
other method, originally developed for the calculation of 
surface elec tronic structure, that is described in full de- 
tail in Refs. Il8ll9l 



In sections III Al and III F3I it was shown that the basic 
quantities that need to be calculated in order to describe 
the scanning probe experiments are transmission coeffi- 
cients from one lead to another. In the next sections, it 
will be shown that these transmission coefficients, and as 
a matter of fact all other quantities we wish to calculate 
(density of states, current density distribution,...) can be 
expressed in terms of the Green's function matrices G^, 
G°„, G„ N and G° Nn (for all possible values of n). These 



M x M matrices G° ; , arc defined as: 



where 



(m\G%\m!) = (m,l\G°\m',l'), 



G° 



(7) 



(8) 



is the Green's function of the device without the influence 
of the tip. In this expression, H d is the Hamiltonian of 
the central device disconnected from the leads, while Sl 
and Sfl are self-energies of the left respectively right lead. 
The Green's functions G^, are thus submatrices of G° 
connecting points from column I' to those of column I of 
the lattice (see also Fig.[Q. 

The method that is used to calculate these Green's 
functions will be presented in Appendix ^ the main 
result being that they can be obtained with a number 



B. Local Scatterer Method 

In our calculations, the scattering potential created 
by the STM tip is modelled by a delta function lo- 
cated on a site (m, n), so it adds a repulsive contribu- 
tion y tlp = v\m,n)(m,n\ to the on-site energies in the 
Hamiltonian (0J. Like explained in section III Al one has 
to calculate a conductance difference: 



Ag(m, n) 



2c 2 



[T° -T(m,n)}, 



(9) 



for every tip position (m, n) in order to obtain the spatial 
map of electron flow in the sample. The transmittancc 
T(m, n) between the left and right lead can be calculated 
b}± 



T(m, n) = Tr [T R G m T L GjyJ, 



(10) 



with ^l(r) related to the retarded self-energies ^l(r) 
the left (right) lead as: 



of 



(11) 



and Gjvi is the Green's function between columns 1 and 
N. The matrix G^vi includes the effect of the STM tip, 
and is therefore different from G^ 1; mentioned in the 
introduction. 
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In Ref.0, one uses the standard recursive Green's func- 
tion method for calculating Gmi'- first, the system is di- 
vided into its separate columns (by making the hopping 
matrices between them zero), and the repulsive potential 
V tlp is added on a certain site. Subsequently the columns 
are attached one by one making use of Dyson's equation. 
For this attachment procedure, N — 1 inversions of an 
M x M matrix are needed to obtain Gni for a single tip 
position. Since this procedure has to be started over and 
over again for every single position of the tip, the total 
number of inversions needed to calculate Gmi for all po- 
sitions scales as MN 2 (in the limit of large N), which 
translates into a total computational effort scaling like 



Since V tlp has only one non-zero element, the inversion 
(1 — G^V^Jf) -1 will boil down to the inversion of a 
scalar. This means that no extra matrix inversions are 
needed to find Gni for an arbitrary lattice position of 
the tip, once we have the G Nn , G°„ and G° x for all n\ 
In Appendix^ we will show a way of calculating these 
functions with a number of matrix inversions that scales 
linear with N, so that the computational effort for cal- 
culating Gni for all tip locations scales as M 3 N. 

At first sight, it might seem that this efficiency is de- 
creased for the complete calculation because one has to 
evaluate the trace in Eq. (|10|) for every tip position, which 
involves products of M x M matrices. The computational 
effort for doing such a product scales like M 3 , and since 
there are MN lattice sites, the total effort would scale as 
M 4 N. However, we have a better way of evaluating this 
trace, scaling as M 3 N, so that we do not loose our effi- 
ciency. Technical details will be described in AppendixlBl 

It now is clear that our method, scaling like M 3 N, 
is more efficient than the standard recursive technique, 
which scales like M 4 N 2 for the same problem. 



C. Voltage Probe Method 

In this case, the STM tip will be modelled by a one- 
dimensional semi-infinite lead, attached to the central 
device at position (m, n) (the lead can be thought to 
extend in a direction perpendicular to the 2D sample). 
The voltage on the tip can be written as a function of 
the transmission coefficients T t i Pi L and T"tip,_R between the 
STM tip and the left and right leads (Eq. ©). These 
transmittances, with the tip positioned over site (m, n), 
can be expressed as: 

Ttip.L = TrTtipG^rLG^], (14a) 
Ttip.jj = Tr[r tip G nN T R Gl N ], (14b) 



M 4 N 2 (computational effort to invert an M x M matrix 
scales as M 3 ). 

However, supposed that we have access to the Green's 
functions G^, G° Nn and G°„ (see Appendix lAl . we are 
able to do things more efficient, by including the effect 
of the tip with Dyson's equation: 

G = G° + G°VG, (12) 

where V is the potential introduced by the tip, V = V tip . 
Projecting l|12l) between columns 1 and N, one obtains (it 
is assumed that the tip is located at lattice site (m,n)): 



(13) 

I 

where r t i P = i (Sy p — S| ip ). Since the lead modelling the 
tip is one-dimensional, one does not have to take care 
of magnetic field effects in this lead, and the self-energy 
S t ip of the tip is known analytically 1 : 

E tip = -te^ Acos ^- E ^ 2t ^. (15) 

The standard recursive Green's function technique 
only gives access to Gni, so it cannot be used to obtain 
the Green's functions G„i and G n N in expression i|14fl . 
But once again, we can use Dyson's equation (fL^)) to re- 
late these Green's functions to Green's functions of the 
device without the tip (calculated in Appendix 1X1) : 

G nl = (1 - G^V^)- 1 G° nl , (16a) 
GnN = (1 - G°„ V^y 1 G nN , (16b) 

where the tip influence leads to an imaginary potential 
V tip = Stip \mn){mn\. 

Again, the inversion (1 — G°„ V^JP)^ 1 will reduce to 
the inversion of a scalar, because V tlp has only one non- 
zero element. In this case, calculation of the traces in 
Eqs. (|14J) is also not computationally expensive since Ttip 
has only one non-zero element. Therefore the total com- 
putational effort scales like M 3 N, needed for the calcu- 
lation of the Green's functions without the influence of 
the tip (see Appendix El . 



D. LDOS - Electron Density 

In order to compare the visualizations of electron flow 
obtained by the two imaging methods above with the 
theory of electron transport, it can be useful to calcu- 
late the local density of states (LDOS) and the electron 
density distribution in the sample, of course without any 
STM tip over the sample. For the LDOS, the standard 



I 

Gni = G Nl + G° Nn (1 — G°„ V^n) 1 G 
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expression isi: 

p(E,m,n) = --Im(m\G° nn (E)\m). (17) 

Having calculated the linear density of states, it is also 
easy to obtain the electron density in the sample by in- 
tegrating over energy. 

E. Current Density Distribution 

When a magnetic field is present, persistent currents 
are flowing through the device, even in the absence of an 



(m,n — 1) — >(m,ra) 



jeq 

(m,n) — >(m+l,n) 



Here, Ep is the Fermi energy of the device, and e is the 
negative electronic charge. We have also introduced: 

S n ft = Vn,n-1 Gn-l,n-l V n -l,rn (19) 

where V n>n -\ describes the hopping between columns n — 
1 and n, and the Green's functions are defined in 

Fig- ID of Appendix |X| 

It is clear, by taking the trace over the row indices 
to, that the total current flowing through every single 
column is equal to zero, so as expected in equilibrium 
there will be no net current through the leads. Also, 
when no magnetic field is present, all Green's functions 



(m,n— 1) — >(m,n) 

A(E) = 

(m,n)— v(m+l,7i) 

C(E) = 



where V = Vl — Vr is the potential difference between 
the leads, and: 

^rj ft = K*,n— 1 ^n— i,n— 1 ^n— l,n> (21a) 
fright _ V n , n +iG%+ 1 ^ n+1 V n +l,n. > (21b) 

T? ht = i[E^ ht -(E^ ht )t]. (21c) 



applied bias. In a recent paper by Cresti et all**, an ex- 
pression for this equilibrium current is derived from the 
Keldysh formalism. Adapted to our notation, the expres- 
sion for the particle current at temperature T = flowing 
from one node to a neighboring node reads (remember 
that to labels the rows of the lattice, n the columns): 



(18a) 
(18b) 

I 

m the Eqs. (HHJ are symmetric so that the equilibrium 
current density in this case vanishes like it should. 

In the non-equilibrium situation, there are two contri- 
butions to the current density. A magnetic field gives 
rise to persistent currents, while the applied bias leads 
to a transport current. Since the persistent currents are 
anti-symmetric with respect to the direction of the field, 
we can define a transport current as the symmetric part 
of the total current density distribution. This transport 
current is gauge invariant, and corresponds to a physi- 
cally relevant (and measurable) quantity. In the linear 
response regime, it is given by (see also Ref . I20T) : 



(20a) 



(20b) 

I 

The Green's functions Gnn in these expressions are 
defined in Fig. of Appendix 1X1 



I 

= 2Re(m| Kn^-^G° nn ] \m), 

= f fj ^ 2Re (to + 1| [C„ (G° nn - G° + )] \m). 



— — {A(E)\ +B + A{E)\_ B } , with 
27m 

2Im(TO| [G°„r^ ht G°t (E left )t] )m) 

{C(E)\ +B + C(E)\_ B }, with 
2Im(TO + l|[C„G°„r^ ht G°t] \m) 
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V. APPLICATIONS 

In the previous sections, we have obtained quite a lot 
of quantities that can give relevant information about the 
flow of electrons in a two-dimensional electron gas. To 
show the power of our method, we discuss two physical 
systems in this section. 



A. Single Quantum Point Contact 

The first system that we consider is the one that is 
used in the experiment of Topinka et ali It consists of a 
2DEG with a quantum point contact in the middle. The 
point contact is modelled by a potential of the form: 

J^exp-^y. (22) 

We have taken a lattice of 1001 by 351 sites, with a lat- 
tice parameter of a = 6.2 nm, which corresponds to a 
hopping parameter t = 14.5 meV (the effective mass of 
the electron is taken to be that for electrons in GaAs: 
rn* = 0.068 m). The parameters of the potential are cho- 
sen as: W = 0.56 1 and £ = 10 a. The Fermi energy is 
put equal to Ep = Lit = 16meV (corresponding to a 
wavelength Xp = 6 a), which is on the first conductance 
plateau of the point contact. For the calculations where 
the tip functions as a local scatterer, the strength of the 
scattering potential is chosen to be v = St. Disorder in 
the system is modelled by a plane of impurities above the 
2DEG, where the repulsive potential from a single impu- 
rity is taken to vary with distance r as l/r 3 , which is 
characteristic for the screened potential in a 2DEG from 
a point charge^. The concentration of impurities is fixed 
at 1% of the total number of lattice sites. The impu- 
rity lattice is located at a distance 6 a above the 2DEG. 
Within the Born approximation, the mean free path of 
the potential we use is estimated to have a value of 4-10 3 a 
which is much longer than the system size so that we are 
in the ballistic regime. The mobility corresponding to 
these parameters is of the order /j, « 2 • 10 6 cm 2 /Vs. 

In Fig. [2^, the calculated current density is seen to 
exhibit a branching behavior that was also apparent in 
the experiment in Ref. 0. The branched flow is also 
present in the map of the conductance difference versus 
the tip position, when the tip is used as a local scatterer 
(Fig. [2J}). It is clear that there is a direct correspon- 
dence between the current density calculation and the 
conductance map. Therefore one can conclude that the 
experiment in Refs. |6j|3 really probes the current distri- 
bution in the sample. Also visible in the conductance 
difference map are the interference fringes spaced by half 
the Fermi wavelength (see the inset in Fig. 12b), which 
were explained as resulting from scattering between the 
point contact and the STM tip&. The map of the lo- 
cal chemical potential, as measured by the STM voltage 
probe (Fig. |2t) in this case gives similar information as 
the previous plots: on the left the current flow appears 
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FIG. 2: Maps of electron flow through a quantum point 
contact. Lengths on the axes are given in units of the lat- 
tice parameter a. Current density distribution (a) in units of 
2e 2 V/(ha), STM conductance map (b) (conductance is mea- 
sured in units of 2e 2 /h) and STM volt probe map (c). The 
voltage on the left lead Vl = 0, and voltage is measured in 
units of Vr. The Fermi energy contour is depicted as a dot- 
ted white line. There is a clear correspondence between the 
different images. 

as regions with increased voltage compared to that of the 
left lead (Vl is put equal to zero). This corresponds to 
a decreased chemical potential due to a deficit of elec- 
trons resulting from the non-equilibrium transport pro- 
cess. On the right, the current flow appears as regions 
with a decreased voltage compared to the right lead. This 
corresponds to an increased chemical potential (excess 
electrons due to the transport process). 

Small oscillations of the chemical potential with a 
wavelength on the order of Xf/2 are apparent in Fig.|2b 
They result from interference between paths which 
emerge from the leads and directly enter the probe, and 
paths which first pass the probe, are reflected from the 
QPC and only then enter the probe. This effect was al- 
ready described in Ref. |U the voltage measurement we 
make is phase sensitive. 

In Fig.|Hl the quantities are calculated for the same sys- 
tem as before, but now a magnetic field is included. The 
field is characterized by a magnetic length of Ib = 28 a, 
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FIG. 3: Maps of electron flow through a quantum point con- 
tact with a moderate magnetic field. Lengths on the axes are 
given in units of the lattice parameter a. Transport current 
density distribution (a) in units of 2e 2 V/(ha), STM conduc- 
tance map (b) (conductance is measured in units of 2e 2 /h) 
and STM volt probe map (c). The voltage on the left lead 
Vl = 0, and voltage is measured in units of Vr. The Fermi 
energy contour is depicted as a dotted white line. 



and a cyclotron radius r c — 835 a. From the non- 
equilibrium transport current density plot (which is sym- 
metric in the magnetic field by its definition in Eas. I2U|) . 
it is clear that the branches of electron flow start bending. 
The radius of curvature has the same order of magnitude 
as the cyclotron radius, so we are seeing here the on- 
set of the skipping orbit movement of the electrons. The 
branches are reflected on the upper and lower edges of the 
sample, a proof that one is still in the ballistic regime. 

The conductance difference map (Fig. \3]p) is quite un- 
clear. This can be explained as resulting from the re- 
duction of backscattering in the presence of a magnetic 
fieldii. But nevertheless the tendency of the branches to 
curve can be observed. 

In this regime, the voltage probe method gives bet- 
ter results. The curved branches are clearly visible in 
Fig. Please keep in mind that the voltage method is 
not symmetric under reversal of the field, which results 
in the asymmetry of the voltage map. This asymmetry 
will be explained in more detail with the help of Fig. 0Ji. 
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FIG. 4: Maps of electron flow through a quantum point 
contact with a high magnetic field. Lengths on the axes are 
given in units of the lattice parameter a. Transport current 
density distribution (a) in units of 2e 2 V/(ha), STM conduc- 
tance map (b) (conductance is measured in units of 2e 2 /h) 
and STM volt probe map (c). The voltage on the left lead 
Vl = 0, and voltage is measured in units of Vr. The Fermi 
energy contour is depicted as a dotted white line. 



In Fig. 01 results are shown for a high magnetic field 
(magnetic length Ib — 4.8 a, cyclotron radius rc = 24 a). 
In a plot of the transport current density the electrons 
are seen to describe skipping orbits along the edges of 
the sample; we are in the quantum Hall regime. In this 
regime the original tip scattering method fails because 
of a lack of backscattering 11 . This is visible in Fig. 0Jd: 
only in the middle of the quantum point contact there is 
some conductance decrease because in this region waves 
travelling in opposite directions are "forced" to overlap. 
The map of the local chemical potential (Fig. gives 
better results: the skipping orbits are clearly visible. 

The asymmetry of this plot was already pointed out 
above, and can be understood as follows: the voltage on 
the right lead is chosen to be higher than that on the 
left lead, so electrons are flowing from the left to the 
right. The magnetic field for this plot is pointing out of 
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FIG. 5: Interference between two crossing gaussian beams. 
The wavelength of the interference pattern depends on the 
angle between the beams. Fringes are more closely spaced for 
an angle of 7r/2 (a) than for an angle 7r/4 (b). 




-2 2 -2 2 5 S 2 

k x k x k x 

FIG. 6: Fourier transforms of the flow maps without a 
magnetic field: local scatterer method (a), voltage probe 
method (b), and current density distribution (c). Wavevec- 
tors are in units of 1/a, the color scale has arbitrary units. 



the plane of the paper, so electrons emerging from the left 
lead flow along the upper left edge of the sample, and this 
edge is in equilibrium with the left lead (no skipping or- 
bits are seen on this edge, only a uniform potential distri- 
bution) . Some of these electrons are transmitted through 
the point contact, which results in a higher chemical po- 
tential (so lower voltage!) than at the upper right 
edge of the sample. The electrons reflected from the con- 
tact (continuing their path on the lower edge) give rise 
to a chemical potential that is lower (= voltage that is 
higher) than fi^ (Vl) on the lower left edge. 

B. Two Quantum Point Contacts 

Looking back at the Figs. [21 and 03 another interest- 
ing interference effect is taking place, which has not been 
observed in the experiment. When the branches of elec- 
tron flow hit the upper and lower border of the sample 
(in the regions from to 200 a and from 800 a to 1000a 
on both figures), there are clear interference fringes vis- 
ible, perpendicular to the border. The wavelength of 
these fringes is larger than that of the fringes observed 
in the scatterer experiment (which resulted from back- 
and forth scattering between the tip and the QPC). This 
interference pattern can be explained as a crossing of two 
or more coherent electron beams (branches). In Fig. EJ 
a simulation is shown where the current density due to 
two crossing gaussian beams with wavevectors ki and 
k2 is calculated. A clear interference pattern is visible, 
extending in the direction ki — k2- From comparison be- 
tween Figs. [5^. andEJ}, it is clear that the wavelength of 
the fringes depends on the angle between the two beams. 
It can be shown that this wavelength is given by: 



with ||ki|| = 1 1 Ids 1 1 = 27r/Ai?, and 9 the angle between the 
beams. 

The different periodicities in the different flow maps 
can be made more visible by doing a Fourier transfor- 
mation. In Fig. this is done for columns 1 to 400 of 
the flow maps in Fig. [21 without a magnetic field. In both 



maps where an STM is used, a circle with radius m 2n/3a 
centered on (k x , k y ) = (0, 0) can be seen. This circle cor- 
responds to the interference pattern resulting from a su- 
perposition of paths between the STM tip and the quan- 
tum point contact, which creates the fringes spaced at 
half the Fermi wavelength. In the current density distri- 
bution (Fig. El;), this circle is of course absent. 

Another feature in Fig.^is the presence of two smaller 
circles centered on the X-axis. They result from the in- 
terference effect between crossing beams explained above. 
Using that the interference pattern of two coherent beams 
is directed along ki — k 2 , together with Eq. [23 these cir- 
cles can indeed be reproduced by having interference be- 
tween a main beam directed along the X-axis and others 
crossing it. While this effect has nothing to do with scat- 
tering off the STM tip, these circles are visible in Fourier 
transforms of all flow maps, including that of the current 
density distribution. 

Now, at crossings between two incoherent beams, an 
interference pattern like in Fig. [SI will not occur. If one 
looks at the plot of the transport current with a mag- 
netic field (Fig.[3J), one can see that some branches bend 
upwards, while other bend downwards under influence of 
the magnetic field. This can be interpreted as follows. In 
the device, the chemical potential will be somewhere be- 
tween that of the left and the right lead. If one assumes 
that the chemical potential on the left lead is larger than 
that on the right, we have an excess of electrons flowing 
from left to right. On the other hand we have a deficit 
of electrons ( "holes" ) flowing the other way. These elec- 
trons and holes bend in opposite ways under influence 
of a magnetic field, because they fill different scattering 
states. It is clear then that electrons and holes, and thus 
branches curving upwards and downwards, are emerging 
from different reservoirs and are thus phase incoherent. 
As a result, one does not expect to see interference be- 
tween beams with different chirality. 

In order to test this statement, we have made some cal- 
culations for a system of two QPC's placed above each 
other. The potential of both QPC's has the same form as 
in Eq.[221 and all parameters are chosen as in Fig. [3 The 
results are shown in Fig. Under influence of the mag- 
netic field, the branches are curved. One can clearly see 
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FIG. 7: Maps of electron flow through a double quantum 
point contact with a moderate magnetic field. Lengths on the 
axes are given in units of the lattice parameter a. Transport 
current density distribution (a) in units of 2e 2 V/(ha), STM 
conductance map (b) (conductance is measured in units of 
2e 2 /h) and STM volt probe map (c). The voltage on the left 
lead Vl = 0, and voltage is measured in units of Vr. The volt- 
age map is symmetrized with respect to the direction of the 
magnetic field. In the conductance map, interference fringes 
with wavelength Af/2 resulting from scattering between STM 
tip and QPC are smoothed out. The Fermi energy contour is 
depicted as a dotted white line. 



the interference between beams with the same chirality, 
but there is no interference at the crossing of two beams 
with opposite chirality, like we expected. This distinc- 
tion becomes clear when comparing the crossings that 
are encircled in Fig. To make things more clear, we 
have smoothed out the interference fringes in Fig. 0j) re- 
sulting from scattering between tip and QPC which had 
a wavelength of Ap/2. In Fig. 0i, we symmetrized the 
voltage probe plot with respect to the direction of the 
magnetic field. Also in this plot, the behavior for coher- 
ent beams crossing is different from that for incoherent 
branches. This proves that the effect could be studied 
experimentally. 



VI. CONCLUSION 



We have established an efficient tight-binding method 
to numerically calculate spatial maps of electron flow as 
obtained in a recent series of scanning probe experiments 
where an STM tip is used as a local scatterer for elec- 
trons 6 . The computational effort of our numerical ap- 
proach scales like M 3 N (in the limit N ^> 1), where 
M is the width of the lattice and N is the length. It 
is in this way more efficient than the standard recursive 
Green's function method which scales like M 4 N 2 for the 
same problem. 

We have also shown expressions for the local density of 
states, the electron density and the current density dis- 
tribution. These quantities cannot be calculated within 
the standard recursive approach, but within our scheme 
they can be expressed in terms of the same Green's func- 
tions already known from the numerical simulation of the 
scanning probe experiment. The computational effort for 
these quantities also scales as M 3 N. 

When a magnetic field is applied, backscattering of 
electrons will be strongly reduced because of the pres- 
ence of edge states so that the original STM method does 
not give the desired results. Therefore, a probe method 
was proposed where the tip is used to measure the local 
chemical potential. For this problem again, the numeri- 
cal effort scales like M N. The image one obtains from 
such a method is not always directly related to the cur- 
rent pattern in the sample, but one can expect to obtain 
relevant information about transport in the sample. 

The power of the method was proven in example cal- 
culations, where a tight-binding lattice with more than 
350000 sites has been used. Moreover, by direct compar- 
ison between a numerical simulation of the experiment 
and a calculation of the exact current density distribu- 
tion, it became clear that the original scanning probe 
technique of Topinka et al& is really imaging current 
flow. Furthermore, a new interference phenomenon has 
been observed which resulted from the crossing of phase 
coherent branches, and a new setup with two QPC's has 
been discussed to distinguish between crossings of coher- 
ent branches and incoherent ones. This distinction is vis- 
ible both when the tip is used as a scatterer, and when 
it is used as a voltage probe, so that an experimental 
investigation of the effect should be possible. 

It should be clear that the method proposed is very 
general, and the information obtained by the different 
imaging tools very broad, so that it can be used to study 
electron flow in a variety of systems ranging from e.g. 
the quantum Hall effect to quantum chaos in electron bil- 
liards^. Moreover, including the spin degrees of freedom 
proves to be rather easy; every matrix element should be 
replaced by a spinor. As such, an even broader range 
of phenomena could be studied, ultimately also those in- 
cluding spin-orbit coupling (e.g. the spin Hall effectS^Si). 
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FIG. 8: Green's function method. After the hopping matrices between columns are set to zero to separate the device in 
vertical slices, and their Green's functions Gf" 1 are calculated (a), the slices are added one by one to obtain (b). Then, 

the isolated slices are added to each other again but now starting from the right (c). Finally, by combining sections pairwisc, 
one can obtain the G% n , G° n and (d). In the figure, the first and last column of the device are colored grey to denote the 
influence of the self-energies of the leads. 



APPENDIX A: CALCULATION OF GREEN'S 
FUNCTIONS RELEVANT TO OUR PROBLEM 

In the paper, it became clear that indeed all quantities 
we need can be expressed in terms of a small subset of 
Green's functions G° x , G°„, G^ N and G° Nn (see Fig.[Q. 
In this appendix, we will treat in detail how to obtain 
these functions with a numerical effort that scales like 
M 3 N. 

The first step is to divide the central device into its sep- 
arate columns, and to put the hopping matrices between 
them equal to zero so that the columns become isolated, 
as depicted in Fig [Ht- Next, one calculates the Green's 
function Gf" 1 for every isolated column i = 1,2,..., N, 
which amounts to doing a single M x M matrix inver- 
sion. This first step thus needs a total of N inversions. 

The next step is to attach the isolated columns one 
by one to each other by including the hopping matrices 
between them using Dyson's equation (this is like in the 
standard recursive technique, see e.g. Ref.0). As such, 
one calculates the Green's functions G°' x and G%£ ( see 
Fig. [SJd) . Please note that since not all columns are at- 
tached to each other in G°' x and G°',f (only those to the 
left of column n), these Green's functions are not equal 
to the Green's functions G^ and G°„ we are trying to 
obtain (an exception of course is G°^ !1 which in fact is 
equal to G 1 ^). The superscript L is added to make this 
distinction clear. 

Attachment of a single isolated column costs one in- 
version of an M x M matrix, so a total of N — 1 matrix 
inversions are necessary to calculate the G^ and G%£ 



for all n, starting from the Gf° l . 

For the third step, we start over from the isolated 
Green's functions already calculated in step 1, and glue 
them together like we did in the previous step on the basis 
of a Dyson's equation, but now beginning from the right, 
like depicted in Fig. [SJ;. The Green's functions we cal- 
culate with every step are G^J, G^ and G°'^. Again, 
a single matrix inversion is needed for the attachment of 
a single column, so a total of N — 1 inversions for the 
completion of the third step. 

The final step is to get the G° Nn , G%, G°„ and G° nl 
we are looking for by attaching the previously calculated 
Green's functions in pairs, like illustrated in Fig. 01. One 
takes a section of connected columns 1 to n (with known 
Green's functions G^ and G^), and attaches it to the 
section of connected columns n + 1 to TV (with Green's 
functions G ^ n+1 and G° n+l n+l ). Please note that a 
single inversion is necessary for each pairwise addition, 
which gives a total of N — 1 inversions for the final step. 

From the previous discussion, it becomes clear that we 
need N + 3{N-1) inversions for calculating G® Nn , G^ N , 
G° n and G^ for all n. So our method scales indeed 
linear with N in the number of inversions. 



APPENDIX B: EVALUATION OF THE TRACE 

For the scatterer method, it is necessary to calculate 
the conductance difference in Eq. For this, it seems 
that we have to evaluate the trace in Eq. ^| for all MN 
tip positions. This trace contains products of M x M 
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matrices, so the numerical effort for this step would scale 
like M 4 N. As such one would loose a factor of M in 
efficiency compared to the rest of the calculation. How- 
ever, there is a better way to evaluate the conductance 
difference in Eq. [5] 

We write (see Eq. iTEfll): 



By substituting Eq. i|Bl|) into Eq. I|10(l . one obtains: 



Gni — Gjsn + A, 
with the M x M matrix: 



(Bl) 



A = G Nn C P i 1 ~ G nn ^nn ) G nl- ( B2 ) 

It is important now to make it clear that since has 
only one non-zero element, namely on position (m,m), 
one can write A as a product of a column matrix and a 
row matrix: 

A = [G^Y n ] TO th co j umn r [G nl ] TO th row , (B3) 

with the scalar r given by (v is the magnitude of the 
repulsive tip potential): 

(B4) 



1 - v(m\G° n \m}' 



T(m, n) = T° + 2RcTr [T R AT L (G° N1 ^] + Tr[T R AT L A^}, 

(B5) 

So in order to evaluate the conductance difference 
Ag(m,n), we need to evaluate only the last two terms 
in Eq. I|B5|) . The last term only involves products of an 
M x AI matrix with row or column matrices because of 
the special form of A. The computational effort for this 
term scales thus as M 2 , which corresponds to a total ef- 
fort of M 3 N for all tip locations. Now, since r i (G^ 1 ) t 
is independent of the tip position, it has to be calcu- 
lated only once (with an effort M 3 ). When this matrix 
is known, the trace in the second term also contains only 
products of an M x M matrix with a row or column ma- 
trix, so the computational effort then scales like M 2 , and 
as a result the effort for all tip positions scales like M 3 N 
in the limit of large N. As such, we do not loose our 
efficiency. 
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